library(pacman)
p_load(tidyverse, data.table, edgeR, ggplot2, hrbrthemes, viridis, scales, rmarkdown)
doc_theme <- theme_ipsum(
base_family = "Arial",
caption_margin = 12,
axis_title_size = 12,
axis_col = "black")
curated_names <- fread(
"../../curated_names.tsv")
aba.bed <- fread(
"../../CP046654.1.bed",
col.names = c(
"chromosome",
"left",
"right",
"locus_tag",
"gene_name",
"strand",
"coding",
"completeness"))
aba.counts <- fread(
"../../all_counts_seal.tsv.gz",
header = FALSE,
col.names = c(
"spacer",
"count",
"condition"))
aba.key <- fread(
"../../aba_key.tsv")
aba.design <- fread(
"../../ABA1_experimental_design.tsv",
na.strings = c("NA"))
aba.genome <- aba.key %>%
group_by(locus_tag) %>%
select(
locus_tag,
spacer,
type,
y_pred,
target,
offset) %>%
inner_join(
curated_names,
by = c("locus_tag" = "AB19606")) %>%
rename(AB19606 = locus_tag)
# define the experimental design space to only take into consideration "tubes"
aba.design <- aba.design %>%
filter(experiment == "tube")
# keep only the counts that are in the experimental design space
aba.counts <- aba.counts %>% semi_join(aba.design)
## Joining, by = "condition"
Helper function to convert from dataframe to data.matrix
data.matrix.withnames <- function(x) {
m <- data.matrix(x[,-1], rownames.force = TRUE)
rownames(m) <- as.data.frame(x)[,1]
return(m)
}
Build up components of bottleneck calculation
# https://www.nature.com/articles/nmeth.3253
aba.design <- aba.design %>%
mutate(generations = case_when(
timing == "T0" ~ 0,
timing == "T1" ~ 9,
timing == "T2" ~ 18))
aba.counts.verbose <- aba.counts %>%
inner_join(aba.design) %>%
inner_join(aba.key) %>%
unite("Sample", drug, dose, timing, sep = "_") %>%
unite("Sample", Sample, rep, sep = "_", remove = FALSE) %>%
select(type, spacer, count, Sample, Sample, generations) %>%
group_by(type) %>%
nest %>%
mutate(
data = map(
data, ~pivot_wider(
.x,
id_cols = spacer,
names_from = Sample,
values_from = count,
values_fill = 0)),
data_matrix = map(
data, ~data.matrix.withnames(.x)))
aba.group <- aba.design %>%
as_tibble %>%
unite("Sample", drug, dose, timing, sep = "_") %>%
pull(Sample) %>%
factor
aba.permut <- model.matrix( ~ 0 + aba.group)
colnames(aba.permut) <- levels(aba.group)
aba.edgeR <- aba.counts.verbose %>%
mutate(
y = map(
data_matrix,
~DGEList(counts = ., group = aba.group, genes = row.names(.))),
keep = map(
y,
~filterByExpr(y = ., design = aba.permut, group = aba.group)),
y = map2(
y,
keep,
~.x[.y, , keep.lib.sizes = FALSE]),
y = map(
y,
~calcNormFactors(.)),
y = map(
y,
~estimateDisp(., aba.permut)),
fit = map(
y,
~glmQLFit(., aba.permut, robust = TRUE)))
Under this paradigm, the fitted values should be equal due to model fitting between replicates.
Thus, we can throw away one replicate and just keep rep. 1 to exhibit the fitted CPM.
aba.cpm <- aba.edgeR %>%
mutate(
cpm.fit = map2(
fit,
y,
#~cpm(y = .x$fitted.values)),
~cpm(y = .x$fitted.values, lib.size = exp(getOffset(.y)))),
cpm.fit = map(
cpm.fit,
~as_tibble(., rownames = "spacer")),
cpm.fit = map(
cpm.fit,
~pivot_longer(.x, !spacer, names_to = "Sample", values_to = "Counts per Million (Fitted)")),
cpm.fit = map(
cpm.fit,
~filter(.x, Sample %in% c(
"None_0_T0_1",
"None_0_T1_1",
"None_0_T2_1"))),
cpm.fit = map(cpm.fit, ~.x %>% mutate(Sample = gsub("_[0-9]$", "", Sample))),
cpm.fit = map(cpm.fit, ~.x %>% select(Sample, spacer, `Counts per Million (Fitted)`)),
cpm.fit = map(cpm.fit, ~.x %>% mutate(Sample = case_when(
Sample == "None_0_T0" ~ "Uninduced",
Sample == "None_0_T1" ~ "Induced for 18 hours",
Sample == "None_0_T2" ~ "Induced for 36 hours"))),
cpm.fit = map(cpm.fit, ~.x %>% mutate(Sample = factor(
Sample,
levels = c(
"Uninduced",
"Induced for 18 hours",
"Induced for 36 hours")))),
Type = type,
# Type = case_when(
# type == "control" ~ "Non-targeting Guides",
# type == "mismatch" ~ "Mismatched Guides",
# type == "perfect" ~ "Perfect Guides"),
Type = factor(
Type,
levels = c(
"control",
"mismatch",
"perfect")))
Do stats to normalize by z-score and change to distribution of controls. https://stats.stackexchange.com/questions/46429/transform-data-to-desired-mean-and-standard-deviation
# aba.cpm <- aba.cpm %>%
# mutate(
# summary = map(
# cpm.fit,
# ~.x %>% group_by(Sample) %>%
# summarise(
# mean = mean(`Counts per Million (Fitted)`) %>% round(3),
# med = median(`Counts per Million (Fitted)`) %>% round(3),
# min = min(`Counts per Million (Fitted)`) %>% round(3),
# max = max(`Counts per Million (Fitted)`) %>% round(3),
# sd = sd(`Counts per Million (Fitted)`) %>% round(3))),
# cpm.fit = map2(
# cpm.fit,
# summary,
# ~inner_join(.x, .y) %>% mutate(z = (`Counts per Million (Fitted)` - mean) / sd )))
#
# aba.cpm <- aba.cpm %>%
# mutate(
# cpm.fit = map(
# cpm.fit,
# ~inner_join(
# .x,
# aba.cpm %>% ungroup %>% filter(type == "control") %>% select(summary) %>% unnest(c(summary)),
# by = "Sample")))
#
# aba.cpm <- aba.cpm %>% mutate(cpm.fit = map(cpm.fit, ~.x %>% mutate(`Normalized Counts per Million (Fitted)` = z * sd.y + mean.y)))
aba.cpm %>%
unnest(cpm.fit) %>%
select(Type, Sample, `Counts per Million (Fitted)`) %>%
ggplot(aes(x = `Counts per Million (Fitted)`, fill = `Sample`)) +
geom_density(lwd = 1, alpha = 0.5) +
scale_fill_manual(values = RColorBrewer::brewer.pal(6,"Paired")[c(TRUE, FALSE)]) +
# scale_fill_viridis(
# discrete = TRUE,
# alpha = 0.7,
# direction = 1,
# option = "inferno") +
scale_x_continuous(
trans = scales::pseudo_log_trans(base = 10),
breaks = c(0, 10^(1:6)),
labels = label_number(scale_cut = cut_short_scale())) +
facet_wrap(~ Type) +
doc_theme
## Adding missing grouping variables: `type`
ggsave(
"Induction_Density.svg",
width = 12,
height = 3,
dpi = 600,
units = "in",
device = "svg")
cpm.max <- aba.cpm %>%
ungroup %>%
select(cpm.fit) %>%
unnest(cols = "cpm.fit") %>%
select(`Counts per Million (Fitted)`) %>% max
aba.cpm.t <- aba.cpm %>%
mutate(
plot = map2(
cpm.fit,
Type,
~ggplot(
.x,
aes(x = `Counts per Million (Fitted)`, fill = `Sample`)) +
geom_density() +
ggtitle(.y) +
doc_theme +
scale_fill_viridis(
discrete = TRUE,
alpha = 0.75,
direction = 1,
option = "inferno") +
scale_x_continuous(
trans = scales::pseudo_log_trans(base = 10),
breaks = c(0, 10^(1:6)),
labels = label_number(scale_cut = cut_short_scale()),
limits = c(0, cpm.max))))
print(aba.cpm.t$plot)
## [[1]]
##
## [[2]]
##
## [[3]]
cpm.max <- aba.cpm %>%
ungroup %>%
select(cpm.fit) %>%
unnest(cols = "cpm.fit") %>%
select(`Counts per Million (Fitted)`) %>% max
aba.cpm.t <- aba.cpm %>%
unnest(cols = c(cpm.fit)) %>%
group_by(Sample) %>%
select(`Type`, `Sample`, spacer, `Counts per Million (Fitted)`) %>%
group_by(Sample) %>%
nest %>%
mutate(
plot = map2(
data,
`Sample`,
~ggplot(
.x,
aes(x = `Counts per Million (Fitted)`, fill = `Type`)) +
geom_density() +
ggtitle(.y) +
doc_theme +
scale_fill_viridis(
discrete = TRUE,
alpha = 0.3,
direction = -1,
option = "magma") +
scale_x_continuous(
trans = scales::pseudo_log_trans(base = 10),
breaks = c(0, 10^(1:6)),
labels = label_number(scale_cut = cut_short_scale()),
limits = c(0, cpm.max))))
print(aba.cpm.t$plot)
## [[1]]
##
## [[2]]
##
## [[3]]
Samples that we’re interested in measuring the differences between
conditions.full <- c(
"None_0_T1 - None_0_T0",
"None_0_T2 - None_0_T0",
"Colistin_0.44_T1 - None_0_T0",
"Colistin_0.44_T2 - None_0_T0",
"Rifampicin_0.34_T1 - None_0_T0",
"Rifampicin_0.34_T2 - None_0_T0",
"Imipenem_0.06_T1 - None_0_T0",
"Imipenem_0.06_T2 - None_0_T0",
"Imipenem_0.09_T1 - None_0_T0",
"Meropenem_0.11_T1 - None_0_T0",
"Meropenem_0.11_T2 - None_0_T0",
"Meropenem_0.17_T1 - None_0_T0",
"Meropenem_0.17_T2 - None_0_T0")
conditions.induction <- c(
"Colistin_0.44_T1 - None_0_T1",
"Colistin_0.44_T2 - None_0_T2",
"Rifampicin_0.34_T1 - None_0_T1",
"Rifampicin_0.34_T2 - None_0_T2",
"Imipenem_0.06_T1 - None_0_T1",
"Imipenem_0.06_T2 - None_0_T2",
"Imipenem_0.09_T1 - None_0_T1",
"Imipenem_0.09_T2 - None_0_T2",
"Meropenem_0.11_T1 - None_0_T1",
"Meropenem_0.11_T2 - None_0_T2",
"Meropenem_0.17_T1 - None_0_T1",
"Meropenem_0.17_T2 - None_0_T2")
conditions.all <- c(conditions.full, conditions.induction)
aba.contrasts <- makeContrasts(
contrasts = conditions.all,
levels = aba.permut)
Perform contrasts and find LFC and FDR
aba.edgeR <- aba.edgeR %>%
mutate(
results = map(
fit,
~glmQLFTest(., contrast = aba.contrasts)),
results = map(
results,
~topTags(., n = Inf)),
results = map(
results, ~as_tibble(.) %>% pull(table)),
results = map(
results, ~rename(., spacer = genes)),
results = map(
results,
~pivot_longer(
.,
!c(spacer, F, PValue, FDR, logCPM),
names_to = "contrast",
values_to = "logFC")),
results = map(
results,
~inner_join(
.,
aba.key %>% select(y_pred, original, locus_tag, offset, spacer))),
results = map(
results,
~.x %>% mutate(contrast = gsub("\\.\\.\\.", " - ", contrast))),
results = map(
results,
~.x %>% mutate(contrast = gsub("logFC\\.", "", contrast))),
summary = map(
results,
~.x %>% group_by(contrast) %>%
summarise(
mean = mean(logFC) %>% round(3),
med = median(logFC) %>% round(3),
min = min(logFC) %>% round(3),
max = max(logFC) %>% round(3),
sd = sd(logFC) %>% round(3) )),
summary.paged = map(
summary,
~paged_table(.)))
## Joining, by = "spacer"
## Joining, by = "spacer"
## Joining, by = "spacer"
print(aba.edgeR$type)
## [1] "mismatch" "perfect" "control"
print(aba.edgeR$summary.paged)
## [[1]]
## contrast mean med min max sd
## 1 Colistin_0.44_T1 - None_0_T0 -0.610 0.002 -11.153 2.306 1.423
## 2 Colistin_0.44_T1 - None_0_T1 -0.380 -0.012 -8.492 2.313 0.950
## 3 Colistin_0.44_T2 - None_0_T0 -1.981 0.073 -16.024 4.483 4.045
## 4 Colistin_0.44_T2 - None_0_T2 -0.767 0.000 -13.714 5.034 2.121
## 5 Imipenem_0.06_T1 - None_0_T0 -0.281 0.007 -6.642 1.207 0.797
## 6 Imipenem_0.06_T1 - None_0_T1 -0.050 0.018 -4.735 1.423 0.502
## 7 Imipenem_0.06_T2 - None_0_T0 -1.431 0.057 -15.434 1.821 2.903
## 8 Imipenem_0.06_T2 - None_0_T2 -0.217 0.020 -10.257 3.698 1.170
## 9 Imipenem_0.09_T1 - None_0_T0 -0.282 0.013 -7.589 1.638 0.891
## 10 Imipenem_0.09_T1 - None_0_T1 -0.052 0.020 -5.682 1.872 0.659
## 11 Imipenem_0.09_T2 - None_0_T2 -0.288 0.017 -12.574 4.461 1.478
## 12 Meropenem_0.11_T1 - None_0_T0 -0.242 -0.031 -5.861 1.241 0.640
## 13 Meropenem_0.11_T1 - None_0_T1 -0.012 -0.002 -4.034 1.718 0.324
## 14 Meropenem_0.11_T2 - None_0_T0 -1.506 0.016 -15.434 2.230 3.104
## 15 Meropenem_0.11_T2 - None_0_T2 -0.292 0.001 -9.671 4.238 1.340
## 16 Meropenem_0.17_T1 - None_0_T0 -0.221 -0.020 -6.323 1.701 0.674
## 17 Meropenem_0.17_T1 - None_0_T1 0.009 0.005 -4.360 2.305 0.431
## 18 Meropenem_0.17_T2 - None_0_T0 -1.607 -0.004 -15.654 2.559 3.285
## 19 Meropenem_0.17_T2 - None_0_T2 -0.393 -0.003 -10.836 5.312 1.619
## 20 None_0_T1 - None_0_T0 -0.230 -0.021 -4.960 0.715 0.555
## 21 None_0_T2 - None_0_T0 -1.214 0.040 -14.465 1.169 2.387
## 22 Rifampicin_0.34_T1 - None_0_T0 -0.173 -0.004 -4.519 0.407 0.489
## 23 Rifampicin_0.34_T1 - None_0_T1 0.057 0.025 -2.478 1.618 0.276
## 24 Rifampicin_0.34_T2 - None_0_T0 -1.500 -0.089 -14.416 4.346 3.075
## 25 Rifampicin_0.34_T2 - None_0_T2 -0.286 -0.048 -9.246 5.597 1.599
##
## [[2]]
## contrast mean med min max sd
## 1 Colistin_0.44_T1 - None_0_T0 -0.846 -0.038 -9.907 1.514 2.016
## 2 Colistin_0.44_T1 - None_0_T1 -0.591 -0.105 -6.993 1.593 1.324
## 3 Colistin_0.44_T2 - None_0_T0 -3.172 -1.846 -15.106 5.213 5.039
## 4 Colistin_0.44_T2 - None_0_T2 -1.716 -0.857 -13.531 4.594 2.734
## 5 Imipenem_0.06_T1 - None_0_T0 -0.367 0.062 -6.966 1.594 1.146
## 6 Imipenem_0.06_T1 - None_0_T1 -0.113 0.021 -5.103 1.470 0.729
## 7 Imipenem_0.06_T2 - None_0_T0 -1.856 -0.939 -14.701 3.312 3.715
## 8 Imipenem_0.06_T2 - None_0_T2 -0.400 0.027 -9.852 5.400 1.612
## 9 Imipenem_0.09_T1 - None_0_T0 -0.386 0.049 -7.180 2.036 1.290
## 10 Imipenem_0.09_T1 - None_0_T1 -0.132 0.008 -5.316 1.956 0.958
## 11 Imipenem_0.09_T2 - None_0_T2 -0.542 -0.003 -11.461 3.937 1.987
## 12 Meropenem_0.11_T1 - None_0_T0 -0.271 0.009 -7.674 1.595 0.934
## 13 Meropenem_0.11_T1 - None_0_T1 -0.016 -0.002 -4.006 1.644 0.464
## 14 Meropenem_0.11_T2 - None_0_T0 -2.069 -0.944 -14.701 4.058 3.976
## 15 Meropenem_0.11_T2 - None_0_T2 -0.613 -0.107 -9.852 4.533 1.796
## 16 Meropenem_0.17_T1 - None_0_T0 -0.285 0.022 -6.566 2.098 0.991
## 17 Meropenem_0.17_T1 - None_0_T1 -0.030 -0.025 -4.461 2.172 0.639
## 18 Meropenem_0.17_T2 - None_0_T0 -2.229 -1.188 -14.701 4.982 4.162
## 19 Meropenem_0.17_T2 - None_0_T2 -0.773 -0.245 -10.259 5.915 2.114
## 20 None_0_T1 - None_0_T0 -0.254 0.064 -7.924 1.002 0.830
## 21 None_0_T2 - None_0_T0 -1.456 -0.501 -12.456 2.511 3.071
## 22 Rifampicin_0.34_T1 - None_0_T0 -0.211 0.108 -5.899 0.747 0.788
## 23 Rifampicin_0.34_T1 - None_0_T1 0.044 0.021 -2.744 2.024 0.428
## 24 Rifampicin_0.34_T2 - None_0_T0 -2.009 -1.123 -14.440 5.162 3.828
## 25 Rifampicin_0.34_T2 - None_0_T2 -0.553 -0.194 -8.684 7.175 2.107
##
## [[3]]
## contrast mean med min max sd
## 1 Colistin_0.44_T1 - None_0_T0 -0.007 -0.001 -1.800 1.405 0.156
## 2 Colistin_0.44_T1 - None_0_T1 -0.008 -0.006 -0.991 1.291 0.149
## 3 Colistin_0.44_T2 - None_0_T0 -0.003 0.011 -8.500 3.440 0.395
## 4 Colistin_0.44_T2 - None_0_T2 0.004 0.003 -6.549 3.240 0.345
## 5 Imipenem_0.06_T1 - None_0_T0 -0.002 -0.002 -1.946 0.874 0.127
## 6 Imipenem_0.06_T1 - None_0_T1 -0.002 -0.004 -0.844 0.901 0.100
## 7 Imipenem_0.06_T2 - None_0_T0 -0.006 -0.002 -6.341 2.089 0.297
## 8 Imipenem_0.06_T2 - None_0_T2 0.001 0.002 -2.179 1.869 0.176
## 9 Imipenem_0.09_T1 - None_0_T0 0.000 -0.001 -1.932 1.286 0.147
## 10 Imipenem_0.09_T1 - None_0_T1 0.000 -0.002 -0.935 1.313 0.124
## 11 Imipenem_0.09_T2 - None_0_T2 0.007 0.000 -3.375 2.626 0.270
## 12 Meropenem_0.11_T1 - None_0_T0 0.006 0.004 -1.181 0.548 0.130
## 13 Meropenem_0.11_T1 - None_0_T1 0.005 0.001 -0.517 0.545 0.117
## 14 Meropenem_0.11_T2 - None_0_T0 -0.003 0.010 -5.389 1.799 0.324
## 15 Meropenem_0.11_T2 - None_0_T2 0.004 0.005 -1.483 1.977 0.233
## 16 Meropenem_0.17_T1 - None_0_T0 0.003 -0.002 -1.088 1.112 0.141
## 17 Meropenem_0.17_T1 - None_0_T1 0.002 -0.001 -0.693 0.779 0.125
## 18 Meropenem_0.17_T2 - None_0_T0 0.000 -0.002 -5.659 2.483 0.427
## 19 Meropenem_0.17_T2 - None_0_T2 0.007 0.005 -1.969 2.662 0.346
## 20 None_0_T1 - None_0_T0 0.001 0.002 -1.102 0.677 0.102
## 21 None_0_T2 - None_0_T0 -0.007 0.002 -4.162 1.043 0.214
## 22 Rifampicin_0.34_T1 - None_0_T0 0.002 0.000 -0.552 0.749 0.087
## 23 Rifampicin_0.34_T1 - None_0_T1 0.001 0.002 -0.462 0.958 0.098
## 24 Rifampicin_0.34_T2 - None_0_T0 -0.015 0.017 -9.194 4.165 0.768
## 25 Rifampicin_0.34_T2 - None_0_T2 -0.008 0.026 -9.216 3.588 0.753
Transpose results: Ungroup from type, then regroup results by contrast group
# aba.edgeR.t <- aba.edgeR %>%
# select(type, results) %>%
# unnest(results) %>%
# group_by(contrast) %>%
# nest %>%
# mutate(
# plot = map2(
# data,
# `contrast`,
# ~ggplot(
# .x,
# aes(x = `logFC`, fill = `type`)) +
# geom_density() +
# ggtitle(.y) +
# theme_ipsum() +
# scale_fill_viridis(
# discrete = TRUE,
# alpha = 0.3,
# direction = -1,
# option = "magma")))
#
# print(aba.edgeR.t$plot)
# aba.edgeR.prior <- fread("../../Results/melted_results.tsv.gz") %>%
# filter(condition %in% c("None_0_T1 - None_0_T0", "None_0_T2 - None_0_T0")) %>%
# select(type, condition, LFC, LFC.adj) %>%
# group_by(condition) %>%
# nest %>%
# mutate(
# plot = map2(
# data,
# `condition`,
# ~ggplot(
# .x,
# aes(x = `LFC.adj`, fill = `type`)) +
# geom_density() +
# ggtitle(.y) +
# theme_ipsum() +
# scale_fill_viridis(
# discrete = TRUE,
# alpha = 0.3,
# direction = -1,
# option = "magma") +
# xlim(
# fread("../../Results/melted_results.tsv.gz") %>%
# filter(condition %in% c("None_0_T1 - None_0_T0", "None_0_T2 - None_0_T0")) %>%
# select(type, condition, LFC, LFC.adj) %>% select(LFC.adj) %>% range)))
#
# print(aba.edgeR.prior$plot)
Show replicates and log2 fold change with respect to T0.
aba.counts.stats <- aba.counts.verbose %>%
select(-data_matrix) %>%
unnest(data) %>%
pivot_longer(
!c(type, spacer),
names_to = "condition",
values_to = "count")
aba.counts.stats <- aba.counts.stats %>%
inner_join(
aba.counts.stats %>%
filter(condition == "None_0_T0_1") %>%
ungroup %>%
rename(base_count = count) %>%
select(spacer, base_count)) %>%
mutate(logFC = log2(count/base_count)) %>%
rename(sample = condition)
## Joining, by = "spacer"
aba.counts.stats <- aba.counts.stats %>%
inner_join(
aba.counts.stats %>%
group_by(sample) %>%
filter(type == "control") %>%
summarise(logFC.med.ctrl = median(logFC))) %>%
mutate(logFC.adj = logFC - logFC.med.ctrl) %>%
separate(
sample, c("drug", "dose", "timing", "rep"), "_") %>%
mutate(rep = paste("Rep", rep)) %>%
mutate(dose.with_units = paste0(dose,"ng/μL")) %>%
unite(
"Sample", drug, dose.with_units, timing, remove = F) %>%
group_by(Sample)
## Joining, by = "sample"
logFC.range <- aba.counts.stats %>% filter(!is.infinite(logFC.adj)) %>% pull(logFC.adj) %>% range
aba.counts.stats.reps <- aba.counts.stats %>%
mutate(Type = factor(type)) %>%
arrange(desc(Type)) %>%
nest %>%
mutate(data = map(
data,
~.x %>% pivot_wider(
id_cols = c(Type, spacer),
names_from = rep,
values_from = logFC.adj))) %>%
filter(Sample != "None_0ng/μL_T0") %>%
mutate(
data.finite = map(
data, ~.x %>% filter(is.finite(`Rep 1`) & is.finite(`Rep 2`))),
correlation = map_dbl(data.finite, ~cor(.x$`Rep 1`, .x$`Rep 2`)),
Title = paste0(
"Log2FC: ", gsub("_", " ", Sample),
"(Cor=", round(correlation, 3), ")"))
aba.counts.stats.reps <- aba.counts.stats.reps %>%
mutate(plot = map2(
data,
`Title`,
~ggplot(
.x,
aes(x = `Rep 1`, y = `Rep 2`)) +
geom_abline(intercept = 0, slope = 1, size = 0.5, colour = "red") +
geom_point(aes(color = Type), shape = 20) +
ggtitle(.y) +
doc_theme +
xlim(logFC.range) +
ylim(logFC.range) +
# scale_alpha_manual(values = c(0.5, 0.15, 0.25)) +
scale_color_viridis(
discrete = TRUE,
alpha = c(0.15, 0.25, 0.35),
direction = -1,
option = "viridis") +
theme(
legend.position = "bottom",
legend.key = element_rect(fill = "#ECECEC")) +
guides(colour = guide_legend(override.aes = list(alpha = 1, size = 5)))))
aba.counts.stats.reps %>% pull(plot) %>% print
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
########
aba.counts.stats.reps %>%
filter(Sample %like% "None") %>%
mutate(Title = gsub("Log2FC: None 0ng/μL ", "", Title)) %>%
mutate(Title = gsub("Cor=", "", Title)) %>%
mutate(Title = gsub("\\(", "~(R^2~", Title)) %>%
unnest(data) %>%
ggplot(
aes(x = `Rep 1`, y = `Rep 2`)) +
geom_abline(intercept = 0, slope = 1, size = 0.5, colour = "black") +
geom_point(aes(color = Type), cex = 2) +
facet_wrap(~Title, labeller = label_parsed) +
scale_colour_manual(
values =
c("control" = alpha("grey", 0.5),
"mismatch" = alpha("#1F78B4", 0.5),
"perfect" = alpha("#E31A1C", 0.5))) +
doc_theme +
theme(legend.position = "none")
Violin Plots
aba.counts.violin <- aba.counts.stats %>%
mutate(Log2FC = logFC.adj) %>%
mutate(rep = paste("Rep", rep)) %>%
mutate(dose.with_units = paste0(dose,"ng/μL")) %>%
unite(
"Sample", drug, dose.with_units, timing, remove = F, sep = " ") %>%
filter((Sample %like% "None") & timing != "T0") %>%
mutate(
Sample = gsub("^", "Induction + ", Sample),
Sample = gsub(" \\+ None 0ng/μL", "", Sample),
Sample = gsub("(T[1,2])", "at \\1", Sample),
Sample = factor(Sample, levels = unique(Sample)),
Sample = factor(Sample, levels = rev(levels(Sample))),
Type = factor(type),
Type = factor(Type, levels = rev(levels(Type)))) %>%
ggplot(aes(x = Sample, y = Log2FC, fill = Type)) +
geom_violin(alpha = 1, width = 0.5) +
doc_theme +
scale_fill_manual(values = rev(RColorBrewer::brewer.pal(6,"Paired")[c(TRUE, FALSE)])) +
coord_flip() +
guides(fill = guide_legend(reverse = TRUE)) +
ggtitle("Distribution of Guide Depletion")
print(aba.counts.violin)
## Warning: Removed 46 rows containing non-finite values (stat_ydensity).
melted_results <- fread("../../Results/melted_results.tsv.gz")
interest <- fread("../../interest.tsv", sep = "\t")
melted_results %>%
filter(condition %in% interest$condition) %>%
filter(shift %like% "None") %>%
mutate(
`Fold-change` = 2^LFC.adj,
Type = type, condition = case_when(
condition %like% "T1" ~ "T1",
condition %like% "T2" ~ "T2")) %>%
ggplot(aes(x = `Fold-change`)) +
geom_density(
aes(fill = Type),
alpha = 0.5) +
facet_wrap(facets = c("condition"), ncol = 1) +
doc_theme +
scale_fill_manual(
values =
c("control" = "grey",
"mismatch" = "#1F78B4",
"perfect" = "#E31A1C"))